if !d.name eq 'PS' then begin
    device,xsize=18,ysize=24,yoffset=3
    !p.thick=1 & !x.thick=1 & !y.thick=1
end

@eu_read_bcprofl
@rhoprof
;!p.multi=[0,3,2]
;!p.charsize=2.5
!p.multi=[0,3,1] & !p.charsize=2.0 & !x.margin=[3.75,0] & !y.margin=[3,0.5]

default,nn,5
ntimes = (size(u))[3]
good = (ntimes-nn)+indgen(nn)-1
print,'doing averages over',nn, ' xaver records ...'

ruu = (total((total(u2,1)/m),1)/l)-(total((total(u^2,1)/m),1)/l)
rvv = (total((total(v2,1)/m),1)/l)-(total((total(v^2,1)/m),1)/l)
rww = (total((total(w2,1)/m),1)/l)-(total((total(w^2,1)/m),1)/l)

uu=fltarr(l,m)
vv=fltarr(l,m)
ww=fltarr(l,m)
oox=fltarr(l,m)
ooy=fltarr(l,m)
ooz=fltarr(l,m)
omega=fltarr(l,m)
gradL_r=fltarr(l,m)
gradL_th=fltarr(l,m)
Frs_r=fltarr(l,m)
Frs_th=fltarr(l,m)
divFrs_r=fltarr(l,m)
divFrs_th=fltarr(l,m)
divFrs=fltarr(l,m)

mFrs_r=fltarr(l,m)
mFrs_th=fltarr(l,m)
div_mFrs_r=fltarr(l,m)
div_mFrs_th=fltarr(l,m)
div_mFrs=fltarr(l,m)

quu=fltarr(l,m)
qvv=fltarr(l,m)
qww=fltarr(l,m)
qwv=fltarr(l,m)
qwu=fltarr(l,m)
qvu=fltarr(l,m)
rsinth=fltarr(l,m)
sinth=fltarr(l,m)
rurms2=fltarr(l,m)

urms2t = ruu+rvv+rww
urms2 = (total(urms2t[good]/float((size(good))[1])))

; reynolds stresses
for k=0, l-1 do begin
  for j=0,m-1 do begin
    rsinth[k,j]=rr*r[k]*sin(theta[j])
    sinth[k,j]=sin(theta[j])
    uu[k,j] = total(u[j,k,good])/float((size(good))[1])
    vv[k,j] = total(v[j,k,good])/float((size(good))[1])
    ww[k,j] = total(w[j,k,good])/float((size(good))[1])
    oox[k,j] = total(ox[j,k,good])/float((size(good))[1])
    ooy[k,j] = total(oy[j,k,good])/float((size(good))[1])
    ooz[k,j] = total(oz[j,k,good])/float((size(good))[1])
    quu[k,j]=rho[j,k]*total(u2[j,k,good]- u[j,k,good]^2)/float((size(good))[1])
    qvv[k,j]=rho[j,k]*total(v2[j,k,good]- v[j,k,good]^2)/float((size(good))[1])
    qww[k,j]=rho[j,k]*total(w2[j,k,good]- w[j,k,good]^2)/float((size(good))[1])
    rurms2[k,j]=quu[k,j]+qvv[k,j]+qww[k,j]
    qwv[k,j]=total(rwv[j,k,good]- rho[j,k]*oz[j,k,good]*v[j,k,good])/float((size(good))[1])
    qwu[k,j]=total(rwu[j,k,good]- rho[j,k]*oz[j,k,good]*u[j,k,good])/float((size(good))[1])
    qvu[k,j]=total(rvu[j,k,good]- rho[j,k]*oy[j,k,good]*u[j,k,good])/float((size(good))[1])
;    qwv[k,j]=total(rwv[j,k,good]- rho[j,k]*w[j,k,good]*v[j,k,good])/float((size(good))[1])
;    qwu[k,j]=total(rwu[j,k,good]- rho[j,k]*w[j,k,good]*u[j,k,good])/float((size(good))[1])
;    qvu[k,j]=total(rvu[j,k,good]- rho[j,k]*v[j,k,good]*u[j,k,good])/float((size(good))[1])
  endfor
endfor

nlev=64
thg=where(th GE -85 AND th LT 85)

omega_0=2.*!pi/(28.*24.*3600.)                 ;28 days
omega[*,1:m-2]=(uu[*,1:m-2]/(rsinth[*,1:m-2]))
omega[*,0]=omega[*,1]
omega[*,m-1]=omega[*,m-2]
omega=omega+omega_0

; computing the gyroscopic pumping balance
; LHS
calL = rsinth^2*omega
for j=0,m-1 do begin
 gradL_r[*,j]=xder_curl(calL[*,j],(r*rr))
; gradL_r[*,j]=deriv(r*rr,calL[*,j])
endfor
for k=0,l-1 do begin
 gradL_th[k,*]=xder_curl(calL[k,*],theta)/(r[k]*rr)
; gradL_th[k,*]=deriv(theta,calL[k,*])/(r[k]*rr)
endfor
;lhs=rh2d*(ww*gradL_r + vv*gradL_th)
lhs=rh2d*(ww*gradL_r + vv*gradL_th)
;

 ; RHS
Frs_r=rsinth*qwu
Frs_th=rsinth*qvu
for k=0,l-1 do begin
 divFrs_th[k,*] = xder_curl(sin(theta)*Frs_th[k,*],theta)
; divFrs_th[k,*] = deriv(theta,sin(theta)*Frs_th[k,*])
endfor
for j=1,m-2 do begin
; divFrs_r[*,j]=deriv(r*rr,(r*rr)^2*Frs_r[*,j])/(r*rr)^2
 divFrs_r[*,j]=xder_curl((r*rr)^2*Frs_r[*,j],r*rr)/(r*rr)^2
 divFrs[*,j]=divFrs_r[*,j] + divFrs_th[*,j]/(rsinth[*,j])
endfor
 divFrs[*,0]=divFrs[*,1]
 divFrs[*,m-1]=divFrs[*,m-2]

nlev=64
thg=where(th GE -85. AND th LT 85)
rgood=where(r ge 0.72 and r le 0.94)

xlength=0.28
xinterval=0.05
xpos1=xinterval
xpos2=xinterval+xlength

lev=grange(-max(lhs[3:l-3,thg]),max(lhs[3:l-3,thg]),nlev)
print,'rho <u_m> . grad L', minmax(lhs[3:l-3,thg])
cpos=[xpos1,0.05,xpos2,0.95]
bpos=[0.15,0.36,0.16,0.65]
;bpos=[xpos1,0.055,xpos2,0.07]
contour,smooth(lhs[*,thg],2),x[*,thg],y[*,thg],/fi,nl=64,/iso,title='!7q!8!Ds!N!8<u!Dm!N> !U.!N !MG!X!8L!6',lev=lev,pos=cpos,$
	ystyle=4,xstyle=4
colorbar,range=[min(lev),max(lev)],pos=bpos,/vertical,$
         ytickformat='(E12.1)',yticks=2,ytickv=[min(lev),0.5*(min(lev)+max(lev)),max(lev)],$
         xaxis=0,char=1.5
;colorbar,range=[min(lev),max(lev)],/horizontal,pos=bpos, $
;         xtickformat='(E12.1)',xticks=2,xtickv=[min(lev),0.,max(lev)],xaxis=0,char=1.5
plots,circle(0,0,0.96),thick=4
plots,circle(0,0,0.718),thick=4
oplot,[0,0],[-0.96,-0.718],thick=4,li=0
oplot,[0,0],[0.718,0.96],thick=4,li=0
oplot,[0.718,0.718],[-1.,1.],thick=2,li=1
oplot,[0.9,0.9],[-1.,1.],thick=2,li=1
;
xpos1=xpos2+xinterval
xpos2=xpos1+xlength
cpos=[xpos1,0.05,xpos2,0.95]
bpos=[xpos1,0.055,xpos2,0.07]
;lev=grange(min(divFrs[3:l-3,thg]),max(divFrs[3:l-3,thg]),nlev)
print,'div Frs = ', minmax(divFrs[3:l-3,thg])
contour,-smooth(divFrs[*,thg],2),x[*,thg],y[*,thg],/fi,nl=64,/iso,title='-!MG!X !U.!N !8F!DRS!N!6',lev=lev,pos=cpos,$
	ystyle=4,xstyle=4
;colorbar,range=[min(lev),max(lev)],/horizontal,pos=bpos, $
;         xtickformat='(E12.1)',xticks=2,xtickv=[min(lev),0.,max(lev)],xaxis=0,char=1.5
plots,circle(0,0,0.96),thick=4
plots,circle(0,0,0.718),thick=4
oplot,[0,0],[-0.96,-0.718],thick=4,li=0
oplot,[0,0],[0.718,0.96],thick=4,li=0
oplot,[0.718,0.718],[-1.,1.],thick=2,li=1
oplot,[0.9,0.9],[-1.,1.],thick=2,li=1

xpos1=xpos2+xinterval
xpos2=xpos1+xlength
cpos=[xpos1,0.05,xpos2,0.95]
bpos=[xpos1,0.055,xpos2,0.07]
balance=lhs+divFrs
;lev=grange(min(balance[3:l-3,thg]),max(balance[3:l-3,thg]),nlev)
print,'balance = ', minmax(balance[3:l-3,thg])
contour,smooth(balance[*,thg],2),x[*,thg],y[*,thg],/fi,nl=64,/iso,title='!7d!8<u!Dm!N>.!6grad!8L!6+!6div !8F!DRS!N!6',lev=lev,$
	pos=cpos,ystyle=4,xstyle=4
;colorbar,range=[min(lev),max(lev)],/horizontal,pos=bpos, $
;         xtickformat='(E12.1)',xticks=2,xtickv=[min(lev),0.,max(lev)],xaxis=0,char=1.5
plots,circle(0,0,0.96),thick=4
plots,circle(0,0,0.718),thick=4
oplot,[0,0],[-0.96,-0.718],thick=4,li=0
oplot,[0,0],[0.718,0.96],thick=4,li=0
oplot,[0.718,0.718],[-1.,1.],thick=2,li=1
oplot,[0.9,0.9],[-1.,1.],thick=2,li=1
end

